#include "processImg.h"
#include "pubHead.h"

namespace analysis_img
{
//static double wnKer[9] = { -1.5, 0.5, -1.5,
//0.5, 5, 0.5,
//-1.5, 0.5, -1.5 };
static double wnKer[9] = { -1, -1, -1,  //v3
-1, 9, -1,
-1, -1, -1 };


float binaryThreshold(const cv::Mat& image){
	int up = 100000;
	int window = 10;
	float bl = 100;
	int width = image.cols;
	int height = image.rows;
	std::vector<long> histData(up + 1, 0);
	int hh = 0;
	unsigned int sum = 0;

	for (int i = 0; i < height; i++)
	{
		for (int j = 0; j < width; j++)
		{
			hh = round(image.at<float>(i, j) * bl);
			if (hh > 0 && hh <= up){
				histData[hh]++;
				sum++;
			}
		}
	}

	float temp1 = 0;
	float temp2 = 0;
	unsigned int sumt = 0;
	for (int i = 0; i < up + 1; i++)
	{
		sumt += histData[i];
		if (sumt * 1.0 / sum > 0.5 && temp1 == 0){
			temp1 = i / bl;
		}
		if (sumt * 1.0 / sum > 0.98 && temp2 == 0){
			temp2 = i / bl;
		}
	}

	return static_cast<float>(cv::mean(image)[0] * (0.6 + temp1 / (temp2 - temp1)));//linear or sqrt
}

void imageFilter_AnaImg::impoccess(cv::Mat &image, cv::Mat &bg, cv::Mat &mse, cv::Mat &inoise, float segma, cv::Mat& onoise, unsigned long &n){

	cv::Mat element = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(17, 17));
	morphologyEx(image, bg, cv::MORPH_OPEN, element);
	cv::Mat noise = cv::Mat(image.rows, image.cols, CV_32F);
	cv::subtract(image, bg, image);

	cv::add(bg, 1, bg);
	cv::divide(image, bg, noise);

	float mv = binaryThreshold(noise) * snrth;
	cv::threshold(noise, inoise, mv, 1.0, cv::THRESH_BINARY);

	mse = image.clone();

	cv::Mat m_kernel;
	m_kernel = cv::Mat(3, 3, CV_64F, wnKer);
	cv::GaussianBlur(image, image, cv::Size(3, 3), segma, segma);
	cv::filter2D(image, image, CV_32F, m_kernel);
	element = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(3,3));
	cv::erode(inoise, onoise, element);
	n = cv::sum(onoise)[0];
}

void imageFilter_AnaImg::OtsuThreshold(const cv::Mat &imgGray, int &sth)
{
	int up = 10000;
	int window = 10;
	int down = 1;

	int width = imgGray.cols;
	int height = imgGray.rows;
	std::vector<long> histData(up + 1, 0);
	std::vector<long> histData2(up + 1, 0);
	int hh = 0;
	for (int i = 0; i < height; i++)
	{
		for (int j = 0; j < width; j++)
		{
			if (imgGray.at<float>(i, j) >= 0 && imgGray.at<float>(i, j) <= up){
				hh = round(imgGray.at<float>(i, j));
				histData[hh]++;
			}
		}
	}
	for (int i = 0; i < up + 1 - window; i++){
		int c = 0;
		int jd = i - window > 0 ? i - window : 0;
		for (int j = jd; j <= i + window; j++){
			histData2[i] += histData[j];
			c++;
		}
		histData2[i] /= c;
	}
	int max = -1;
	for (int i = down; i < up + 1; i++)
	{
		if (histData2[i] > max){
			max = histData2[i];
			sth = i;
		}
	}

	float dest = static_cast<float>(max * 0.7);
	float mind = 1e10;
	int end = sth;
	for (int i = end; i < up + 1; i++){
		if (std::abs(histData2[i] - dest) < mind){
			mind = std::abs(histData2[i] - dest);
			sth = i;
		}
	}
	sth = static_cast<int>(end * (1.0 + std::sqrt(end * 1.0 / (sth - end))));//linear or sqrt
}
};